R Training 3

Introduction

Welcome to the third tutorial! In this lesson we are going to look at actual bird songs and learn how identify and seperate syllables. We will also look at some spike data at the end. By this time, you should be familiar with the basics of R, so there will be less syntax explanations. However, if you ever get confused or want to know more about a function, you can search for it using the Help tab ion the bottom right pane.

Note: We will be using the default graphics for some of our plots for conveniance, as spectro() uses default graphics and will not allow mulit-plots with ggplot2 generated graphs. There is a ggspectro() function, but it is complicated and we will not be using it.

Exercise 9: Sperate the syllables in a song

Lets load our data. Remember to load any necessary libraries before you begin.

library(seewave)
library(tuneR)
library(ggplot2)
library(gridExtra)
zf_data = readWave('zfinch.wav')
zfFs = zf_data@samp.rate
tim = seq(1/zfFs, length(zf_data@left)/zfFs, 1/zfFs)

We are going to separate the syllables by amplitude. The problem with the original signal is that it goes both up and down, so it is not possible to directly measure the amplitude. Below I describe two methods,“Strategy 1” and “Strategy 2”, that can be used to make a useful measure of the amplitude.

These two methods produce nearly the same result. The first method is to rectify the signal and low-pass filter, and the other is to use a function known as “Hilbert” and low-pass filter.

Both require that the signal be centered at zero. It is easy to do this by subtracting the mean from the original signal.

Lets use the mean() function…

zf_data_center = zf_data@left - mean(zf_data@left)

Strategy 1: Rectify and smooth

First lets take the absolute value using abs()…

rz = abs(zf_data_center)

Now we are going to smooth this signal using a median filter. R does this with the function runmed().

Before we do that lets set our variables. Here sms is the duration, in seconds, of the filter - I have set the default to 20 milliseconds. If sms is too long, then it will smooth out the signal too much. If sms is too short, then you will get spurious syllables. We will multiply sms with the sample rate, Fs. This will give us the number of samples for the duration that sms specifies.

sms = 0.020
zfFs = zf_data@samp.rate

The we apply our filter…

mrz = runmed(rz, Fs*sms)

We don’t need to do the next step, but it makes life a little easier later. Here we multiple by 1000 so that the values are ~ 1 rather than 0.001.

mrz = mrz*1000

Strategy 2: Hilbert and Smooth

This is almost identical to Strategy 1, but uses a nice function called hilbert(). This function is from the package seewave, which you should already have installed!

hz = abs(hilbert(rz, f = Fs))

Lets apply the same filter as above.

mhz = runmed(hz, Fs*sms)
mhz = mhz*1000

Time to plot. First the original signal…

zf_data_plot = data.frame(tim, zf_data@left)
ggplot(zf_data_plot, aes(x = tim, y = zf_data.left)) + 
  geom_line(color = "Blue") +
  scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
  scale_x_continuous(name = "Time")

…and the same signal after filtering.

ggplot() +
  geom_line(aes(x = tim, y = mhz), color = "Green") +
  geom_line(aes(x = tim, y = mrz), color = "Red") +
  ylab("Time") +
  xlab("Frequency")

Sepperate the Song into Syllables

Animals often produce signals that are composed of multiple parts. One strategy to understand the signal is to divide it into syllables. That is what we will do here - use an amplitude threshold to separate each part into distinct syllables.

We start by replotting the same data (mrz) from the last system, but the y axis will be on a log scale. This is acheived using the parameter scale_y_log10() which plots with the X axis normal and Y axis on a log scale.

ggplot() +
  geom_line(aes(x = tim, y = mrz), color = "Red") +
  scale_y_log10() +
  theme(legend.position = "none")

Now that we’ve plotted our data, need to choose a threshold. We can either look at the plot and guess, or we can click on the plot and the computer will give us the coordinates where we clicked:

First strategy - I guess “900”!

thresh = 900

For clicking, we use the command locator() and specify how many clicks we are going to make, which is 1 in this case:

thresh = locator(1)

If we used locator(), then thresh will have two values, the x-value and the y-value of the click. We only want the y-value, which is in the second position, and hence we only want thresh[2]. We can do the following:

thresh = thresh[2]

Let’s replot with a line at the guessed threshold just to check. You might have to redo this section over and over again with different values for thresh until we get a value that we think will suffice.

And now we plot the line of our threshold :

ggplot() +
  geom_line(aes(x = tim, y = mrz), color = "Red") +
  geom_line(aes(x = tim, y = thresh), color = "Blue") +
  scale_y_log10() +
  theme(legend.position = "none")

NA

Now we use the which() command to get the part of the signal above the threshold. The variable syls will get the index numbers for each value of mrz that is larger than thresh:

syls = which(mrz > thresh)

Now we’re going to do something that will extract the starts and ends of each syllable. This seems complicated, but it is simple.

The first step is to make a list of zeros that is the length of the signal. We can use the rep() command:

zz = rep(0, length((zf_data@left)))

Let’s look at zz:

Replot the original signal and then plot zz.

grid.arrange(ggplot(zf_data_plot, aes(x = tim, y = zf_data.left)) + 
               geom_line(color = "Blue") +
               scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
               scale_x_continuous(name = "Time"),
             ggplot() +
               geom_line(aes(x = tim, y = zz, color = "Red")) +
               theme(legend.position = "none"), 
             nrow = 2)

From above, syls is the list of values that are above the threshold - these are where the syllables are located. We will set zz to 1 for each of these values…

zz[syls] = 1

Now lets replot zz:

grid.arrange(ggplot(zf_data_plot, aes(x = tim, y = zf_data.left)) + 
               geom_line(color = "Blue") +
               scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
               scale_x_continuous(name = "Time"),
             ggplot() +
               geom_line(aes(x = tim, y = zz, color = "Red")) +
               theme(legend.position = "none"), 
             nrow = 2)

Cool - these blocks roughly align with the syllables… but what we need are the start times and end times of the syllables. We will use a trick to get them…

diff takes the difference between adjacent values…

For example sample(2) - sample(1)…

yy = diff(zz)
yy_tim = tim[1:length(tim)-1]

Now let’s plot this. Starts are marked by 1 and ends are marked by -1.

ggplot() +
  geom_line(aes(x = yy_tim, y = yy), color = "Blue")

Now we use the which() command to get the starts and the ends of each syllable. Notice that we use == to indicate when we are ASKING if values are equal, but use = to SET a variable to a value.

starts = which(yy == 1)
ends = which(yy == -1)

Now this is simple! We use a loop - the for command to get each syllable. Before we write our loop we will create some empty lists to hold our outputs.

Note: Loops are generally not recommended in R, but here it is fairly simple and makes ourlives easier.

syllable = c()
timmy = c()
timm = c()
for (i in seq(1,length(starts))) {
  syllable[[i]] = zf_data@left[starts[i]:ends[i]]
  timmy[[i]] = tim[starts[i]:ends[i]]
  timm[[i]] = seq(1/Fs, (1 + ends[i] - starts[i])/Fs, 1/Fs)
}

syllable[1] will have the entire 1st syllable, syllable[2] will have the entire 2nd syllable, etc.

Note: timmy and timm are not necessary, but might be useful. timmy is the real times from the song. timm is the time of the syllable from 0

Just to show what has happened, here is a plot:

par(mfrow = c(2,4))
specplot(syllable[[4]], Fs = zfFs)
specplot(syllable[[5]], Fs = zfFs)

specplot(syllable[[6]], Fs = zfFs)
specplot(syllable[[7]], Fs = zfFs)

plot(timmy[[4]], syllable[[4]], type = "l", col = "Red")
plot(timmy[[5]], syllable[[5]], type = "l", col = "Red")

plot(timmy[[6]], syllable[[6]], type = "l", col = "Red")
plot(timmy[[7]], syllable[[7]], type = "l", col = "Red")

I picked syllables 4,5,6, and 7. You can pick others if you wish.

Now we will get the silent parts between syllables. This can be important because sometimes the time between signals is an independent signal. For example, in frogs the duration between calls determines the “pulse repetition rate”, which can indicate whether the call is a mate attraction signal or an aggressive signal.

The procedure is almost the same as for the syllables as above. What is different is that instead of copying the signal from zfData, we instead make pure silence by putting in a flat signal with an amplitude of 0.

We use a loop for each of the ends (which is the start of each silence).

To do this, make a list of zeros with the length of the interval and get the length of the interval by getting the difference between the end of the silence, which is the ‘start’ of the next syllable, and the start of the silence, which is the ‘end’ of the previous syllable. Complicated!?! Not really!

nop = c()
noptimmy = c()
noptimm = c()
 
for (j in seq(2,length(ends))) {
  
  nop[[j]] = list(rep(127, times = starts[j] - ends[j-1]))
  noptimmy[[j]] = tim[ends[j-1]:starts[j]]
  noptimm[[j]] = seq(1/Fs, (1 + ends[j] - starts[j])/Fs, 1/Fs)
  
}

For review: You can easily get the data for each syllable.

syllable[[1]] is the data for the first syllable, and syllable[[2]] the second.

Let’s plot a syllable.

par(mfrow = c(1,2))
specplot(syllable[[7]], Fs = zfFs, ovlp = 90)
plot(timmy[[7]], syllable[[7]], type = "l", col = "Red")

How many syllables did we find??

length(syllable)

How long (duration) is syllable{7}?

length(syllable[[7]])/Fs

Exercise 10: Make a NEW song with the syllable order reversed

Here we are going to take each syllable, and instead of the normal order of ABCDEFG we will make the order GFEDCBA. Easy!

Here we are making a new variable, revordersong and we are putting the last syllable in there.

revordersong = syllable[length(syllable)]
revordersong = unlist(revordersong)

Now we will cycle though each syllable from the second to last to the first (we already have the last syllable from the above line.

The variables in the brackets will be concatenated. So [ revordersong nop syllable ] will make the previously defined revordersong followed by the silent period, which is then followed by the next syllable.

for (k in seq(length(syllable)-1, 1, -1)) {
  
  revordersong = append(revordersong, unlist(nop[[k + 1]]))
  revordersong = append(revordersong, unlist(syllable[[k]]))
  
}

This is nice to add - a time base for this new signal. We’ve seen this many times by now.

revtim = seq(1/Fs, length(revordersong)/Fs, 1/Fs)

p1 = ggplot() + 
  geom_line(aes(x = tim, y = zf_data@left), color = "Blue") +
  scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
  scale_x_continuous(name = "Time")
p2 = ggplot() +
  geom_line(aes(x = revtim, y = revordersong), color = "Red") +
  scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
  scale_x_continuous(name = "Time")
grid.arrange(p1, p2, nrow = 2)

or we can show the spectra…

the orininal signal:

specplot(zf_data, Fs = zfFs)

reverse signal:

specplot(revordersong, Fs = zfFs)

Save the data to a wave file. You can then use Audacity to listen to it.

savewav(revordersong, Fs, "revordersong.wav")

Exercise 11: Make a NEW song with each syllable reversed but the syllable order normal

This is nearly identical to the previous excerise, but now the order is correct (ABCDEF) but each syllable is reversed.

Take the first syllable and put it in our new variable revsylsong…

revsylsong = unlist(syllable[length(syllable)])

Now reverse it – we saw how to do this in a previous excerise…

revsylsong = revsylsong[seq(to = 1, by = -1, length.out = length(revsylsong))]

Now do the same for the rest of the syllables, adding them to the new variable…

We use the temporary variable ra to do the reversing

for (l in seq(2, length(syllable), 1)) {
  
  ra = unlist(syllable[[l]])
  revsylsong = append(revsylsong, unlist(nop[[l]]))
  revsylsong = append(revsylsong, unlist(ra[seq(to = 1, by = -1, length.out = length(syllable[[l]]))]))
  
}

and don’t forget about the time data…

revsyltim = seq(1/Fs, length(revsylsong)/Fs, 1/Fs)

Now we can make the same plots as last time, including spectrograms:

p3 = ggplot() + 
  geom_line(aes(x = tim, y = zf_data@left), color = "Blue") +
  scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
  scale_x_continuous(name = "Time")
p4 = ggplot() +
  geom_line(aes(x = revsyltim, y = revsylsong), color = "Red") +
  scale_y_continuous(name = "Frequency", breaks = c(27,127,227), limits = c(0,250)) +
  scale_x_continuous(name = "Time")
grid.arrange(p3, p4, nrow = 2)

We can save the file as a .wav using the same process as in the last excercise.

savewav(revsylsong, FS, "revsylsong.wav")

With these tools you can make a wide range of arbitrary stimuli from your recordings of animal signals. These are powerful tolls for sensory electrophysiology.

Practice:Try to redo this exercise with the data in wcs.wav.

Exercise 12: Spikes

In this exercise we will take a recording of a neural activity (spikes!) and turn it into data that we can start to analyze in useful ways.

I have provided two sample files that contain some nice neural activity. They are very different from each other.

Read the first file, make the time data, and plot (easy by now).

a = readWave("spikes1.wav")
aFs = a@samp.rate
atim = seq(1/aFs, length(a)/aFs, 1/aFs)
ggplot() +
  geom_line(aes(x = atim, y = a@left), colour = "Blue")

Lets zoom in on a second of the data…

Note: This can be done several ways. The first way is to us p + scale_x_continuous(limits = c(3,4)). This way will remove all points not plotted. You can use xlim(3,4) as a shorthand version of this same function. The method we use is p + coord_cartesian(xlim = c(3,4)), which PRESERVES data points not plotted. This may be relevant if you are fitting anything to the data.

ggplot() +
  geom_line(aes(x = atim, y = a@left), colour = "Blue") + 
  coord_cartesian(xlim = c(3,4))

Now lets set and amplitude threshold for the spikes. Different recordings have different amplitudesfor spikes. Here I guessed 15000. Lets look at the plot to confirm.

athresh = 15000

We will use the which() function to get all of the data above the threshold.

spikes = which(a@left > athresh)

This is a complicated, but necessary step. Spikes are more than 1 or two milliseconds in duration. Thus, we want a threshold for the interval between spikes - anything less than that threshold is a mistake from the **which()* command.

Our threshold will be 1 millisecond: 0.001 seconds. Since are data is in samples and not in milliseconds, we’ll set our threshold to the sample rate * milliseconds.

diff is the difference between adjacent data points.

So, here we find all of the instances where the number of samples between spikes is greater than our threshold of 1 millisecond samples (which for Fs = 10000 is 10 samples).

spikes1 = spikes[which(diff(spikes) > 0.001 * aFs)]

How many spikes did we find?

length(spikes1)

[1] 223

What is the spike rate?

length(spikes1)/atim[length(atim)]

[1] 44.59911

Now let’s plot. We’ll plot a STAR at each time that a spike occurred. We’ll use the rep() function, except this time we will create a list of ones. This would be OK, but I decided to multiply it by the threshold so that the dots will be at the level of the threshold when we plot it. you will see…

ys = rep(1, length(spikes1)) * athresh

Now plot…

ggplot() +
  geom_line(aes(x = atim, y = a@left), colour = "Blue")

…and add the stars.

ggplot() +
  geom_line(aes(x = atim, y = a@left), colour = "Blue") +
  geom_point(aes(x = atim[spikes1], y = ys), colour = "Red", size = 2, shape = 1)

Zoom in to get a better look!

ggplot() +
  geom_line(aes(x = atim, y = a@left), colour = "Blue") +
  geom_point(aes(x = atim[spikes1], y = ys), colour = "Red", size = 2, shape = 1) +
  coord_cartesian(xlim = c(1,2))

We can go even further…

ggplot() +
  geom_line(aes(x = atim, y = a@left), colour = "Blue") +
  geom_point(aes(x = atim[spikes1], y = ys), colour = "Red", size = 2, shape = 1) +
  coord_cartesian(xlim = c(1,1.25))

Another important measure is the time between spikes. This is very easy to get using the diff() command.

diff() will give use the number of samples between spikes, and then we divide by the sample rate to get the time, in seconds, between spikes.

intervals1 = diff(spikes1)/aFs

Here are some easy measures.

mean(intervals1)

[1] 0.02217613

pmax(intervals1)

  [1] 0.0117 0.0374 0.0065 0.0029 0.0765 0.0307 0.0092 0.0516 0.0251 0.0046 0.0056 0.0524 0.0062
 [14] 0.0702 0.0044 0.0773 0.0031 0.0047 0.0025 0.0197 0.0182 0.0035 0.0033 0.0040 0.0411 0.0046
 [27] 0.0030 0.0287 0.0036 0.0045 0.0579 0.0047 0.0034 0.1463 0.0027 0.0055 0.0037 0.0036 0.0312
 [40] 0.0186 0.0046 0.0046 0.0031 0.0232 0.0035 0.0526 0.0056 0.0028 0.0078 0.0344 0.0058 0.0717
 [53] 0.0028 0.0033 0.0059 0.0543 0.0032 0.0028 0.0612 0.0042 0.0043 0.0025 0.0551 0.0030 0.0031
 [66] 0.0058 0.0241 0.0028 0.0252 0.0024 0.0050 0.0884 0.0061 0.0112 0.0543 0.0089 0.0080 0.0048
 [79] 0.0539 0.0315 0.0044 0.0038 0.0946 0.0349 0.0023 0.0063 0.0039 0.0358 0.0043 0.0273 0.0022
 [92] 0.0070 0.1410 0.0037 0.0351 0.0046 0.0130 0.0044 0.0061 0.0210 0.0026 0.0055 0.0025 0.0369
[105] 0.0300 0.0023 0.0058 0.0073 0.1415 0.0320 0.0695 0.0026 0.0221 0.0045 0.0105 0.0067 0.0045
[118] 0.0682 0.0027 0.0051 0.0027 0.0043 0.0759 0.0051 0.0964 0.0385 0.0035 0.0048 0.0037 0.0582
[131] 0.0025 0.0055 0.0048 0.0815 0.0219 0.0028 0.0050 0.0306 0.0036 0.0027 0.0053 0.0323 0.0074
[144] 0.0042 0.0582 0.0043 0.0550 0.0068 0.0375 0.0027 0.0071 0.0306 0.0277 0.0282 0.0027 0.0031
[157] 0.0684 0.0034 0.0060 0.1088 0.0046 0.0030 0.0038 0.0816 0.0194 0.0022 0.0029 0.0042 0.0920
[170] 0.0031 0.0054 0.0042 0.0029 0.0328 0.0600 0.0656 0.0345 0.0034 0.0050 0.0055 0.0158 0.0041
[183] 0.0719 0.0043 0.0397 0.0036 0.0035 0.0246 0.0076 0.0022 0.0046 0.0032 0.1378 0.0025 0.0462
[196] 0.0043 0.0067 0.0026 0.0055 0.0319 0.0057 0.0942 0.0514 0.0047 0.0057 0.0369 0.0054 0.0412
[209] 0.0159 0.0026 0.0058 0.0031 0.0023 0.0043 0.0567 0.0255 0.0039 0.1530 0.0034 0.0068 0.0033
[222] 0.0082

pmin(intervals1)

  [1] 0.0117 0.0374 0.0065 0.0029 0.0765 0.0307 0.0092 0.0516 0.0251 0.0046 0.0056 0.0524 0.0062
 [14] 0.0702 0.0044 0.0773 0.0031 0.0047 0.0025 0.0197 0.0182 0.0035 0.0033 0.0040 0.0411 0.0046
 [27] 0.0030 0.0287 0.0036 0.0045 0.0579 0.0047 0.0034 0.1463 0.0027 0.0055 0.0037 0.0036 0.0312
 [40] 0.0186 0.0046 0.0046 0.0031 0.0232 0.0035 0.0526 0.0056 0.0028 0.0078 0.0344 0.0058 0.0717
 [53] 0.0028 0.0033 0.0059 0.0543 0.0032 0.0028 0.0612 0.0042 0.0043 0.0025 0.0551 0.0030 0.0031
 [66] 0.0058 0.0241 0.0028 0.0252 0.0024 0.0050 0.0884 0.0061 0.0112 0.0543 0.0089 0.0080 0.0048
 [79] 0.0539 0.0315 0.0044 0.0038 0.0946 0.0349 0.0023 0.0063 0.0039 0.0358 0.0043 0.0273 0.0022
 [92] 0.0070 0.1410 0.0037 0.0351 0.0046 0.0130 0.0044 0.0061 0.0210 0.0026 0.0055 0.0025 0.0369
[105] 0.0300 0.0023 0.0058 0.0073 0.1415 0.0320 0.0695 0.0026 0.0221 0.0045 0.0105 0.0067 0.0045
[118] 0.0682 0.0027 0.0051 0.0027 0.0043 0.0759 0.0051 0.0964 0.0385 0.0035 0.0048 0.0037 0.0582
[131] 0.0025 0.0055 0.0048 0.0815 0.0219 0.0028 0.0050 0.0306 0.0036 0.0027 0.0053 0.0323 0.0074
[144] 0.0042 0.0582 0.0043 0.0550 0.0068 0.0375 0.0027 0.0071 0.0306 0.0277 0.0282 0.0027 0.0031
[157] 0.0684 0.0034 0.0060 0.1088 0.0046 0.0030 0.0038 0.0816 0.0194 0.0022 0.0029 0.0042 0.0920
[170] 0.0031 0.0054 0.0042 0.0029 0.0328 0.0600 0.0656 0.0345 0.0034 0.0050 0.0055 0.0158 0.0041
[183] 0.0719 0.0043 0.0397 0.0036 0.0035 0.0246 0.0076 0.0022 0.0046 0.0032 0.1378 0.0025 0.0462
[196] 0.0043 0.0067 0.0026 0.0055 0.0319 0.0057 0.0942 0.0514 0.0047 0.0057 0.0369 0.0054 0.0412
[209] 0.0159 0.0026 0.0058 0.0031 0.0023 0.0043 0.0567 0.0255 0.0039 0.1530 0.0034 0.0068 0.0033
[222] 0.0082

But looking at the plot we see that the spikes are in “bursts” rather than randomly spread. If they were randomly timed, we would find a random distribution of the intervals between spikes. Lets take a look using a histogram…

ggplot() +
  geom_histogram(aes(x = intervals1), bins = 100)

For these data we see that there are many many intervals below 0.015 (15 milliseconds), and a spread of longer intervals. This plot is a very good example of a “bursty” neuron.

The other data, spikes2.wav, is very very different.

Please do the same analysis for those data as you did for spikes1.wav. Use different variable names like changing spikes1 to spikes2 and interval1 to interval2 so that you can plot a histogram just like we did above. Try plotting both histograms in the same frame.